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Abstract 

Background: The ability to measure and quantify myocardial motion and deformation provides a useful tool to 
assist in the diagnosis, prognosis and management of heart disease. The recent development of magnetic 
resonance imaging methods, such as harmonic phase analysis of tagging and displacement encoding with 
stimulated echoes (DENSE), make detailed non-invasive 3D kinematic analyses of human myocardium possible in 
the clinic and for research purposes. A robust analysis method is required, however. 

Methods: We propose to estimate strain using a polynomial function which produces local models of the 
displacement field obtained with DENSE. Given a specific polynomial order, the model is obtained as the least 
squares fit of the acquired displacement field. These local models are subsequently used to produce estimates of 
the full strain tensor. 

Results: The proposed method is evaluated on a numerical phantom as well as in vivo on a healthy human heart. 
The evaluation showed that the proposed method produced accurate results and showed low sensitivity to noise 
in the numerical phantom. The method was also demonstrated in vivo by assessment of the full strain tensor and 
to resolve transmural strain variations. 

Conclusions: Strain estimation within a 3D myocardial volume based on polynomial functions yields accurate and 
robust results when validated on an analytical model. The polynomial field is capable of resolving the measured 
material positions from the in vivo data, and the obtained in vivo strains values agree with previously reported 
myocardial strains in normal human hearts. 



Background 

The pumping behavior of the heart consists of complex 
sequences that constitute cardiac contraction and 
relaxation. The kinematic behavior of the heart has been 
analyzed extensively in order to understand the mechan- 
isms that impair the contractile function of the heart 
during disease. Until recently, the only method with 
high enough spatial resolution of three-dimensional 
(3D) myocardial displacements to resolve transmural 
behaviors was invasive marker technology [1,2]. How- 
ever, the recent development of magnetic resonance 
imaging (MRI) methods such as harmonic phase 
(HARP) analysis of tagging [3] and displacement 
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encoding with stimulated echoes (DENSE) [4], make 
detailed non-invasive 3D transmural kinematic analyses 
of human myocardium possible for clinical and research 
purposes [5]. 

A previously presented polynomial method for cardiac 
strain quantification from surgically implanted markers 
and beads enables straightforward 3D strain computa- 
tion within the myocardium [6]. For the marker and 
bead data, this method was shown to yield accurate and 
robust results, with errors smaller or comparable to a 
previously presented finite element method tailored for 
the same type of data, and has been applied to bead dis- 
placements for analyses of systolic and diastolic myocar- 
dial strains [7-9]. The method is simple in nature which 
aids to bridge a possible gap in understanding between 
different disciplines and has specifically been shown to 
be well suited for sparse arrays of displacement data [6]. 
The present work applies this polynomial method to 3D 
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MRI displacement data as it: 1) allows for quantification 
of the full 3D strain tensor; 2) resolves transmural strain 
variations within the myocardium, and 3) provide 
robustness to noise. Displacement measurements were 
acquired within a 3D myocardial slab at the left ventri- 
cular (LV) equator. The 3D volume was divided into 
segments where the polynomial method was applied to 
quantify the Lagrangian strain tensor. The method was 
evaluated on a numerical phantom as well as in vivo on 
a healthy human heart. 

Methods 

Strain analysis 

Displacement-encoded MRI acquire displacements u at 
time t n of timeframe n relative to a reference configura- 
tion of the myocardium, at the spatial coordinates s = 
(s x , s y , s z ). Lower-case letters are used to denote the 
coordinates of the deformed myocardium. The spatial 
coordinates S = (S x > S Y > S z ) of the myocardium in the 
reference configuration were derived by subtracting the 
displacements from the spatial coordinates at the cur- 
rent configuration; S(ti) = s - u(s, ti). The upper-case 
letters are used to denote the coordinates of the refer- 
ence configuration. 

X = X(S) is the coordinate of a material point, defined 
as an infinitesimal volume of myocardial tissue, in the 
reference configuration and x = x(s) is the coordinate 
after a deformation of the body. The mapping from 
reference to deformed configuration was modeled by a 
polynomial function g(X) of the coordinate of the mate- 
rial point in the reference configuration [6]. This poly- 
nomial position field gives an estimate x of the 
measured coordinate x and can be of different order in 
the different reference coordinate directions, radial (X R ), 
circumferential (X c ), and longitudinal (X L ), depending 
on the number of material points along each dimension, 
and can thus be described as 



x=g(X) = f R (X R )f c (Xc)f L (X L ) / 



(1) 



where f R , f c and f L are polynomial functions of X R , X c 
and X L , respectively. The Lagrangian strain tensor E is a 
measure of deformation and is connected to the defor- 
mation gradient tensor F via the definition 



E = 0.5(F T F-I), 



(2) 



where I is the identity tensor. The deformation gradi- 
ent tensor is given by differentiation, with respect to 
reference position, of the mapping from reference to 
deformed configuration, F = dx/dX. 

The principal strains E lf E 2 and E 3 are obtained by 
diagonalization of the strain tensor, and their magni- 
tudes are independent of any reference coordinate 



system. Principal strain E 1 represents maximum length- 
ening and E 3 represents maximum shortening. 

In the subsequent analytical and in vivo evaluation, 
the 3D volumes were divided into 80 overlapping seg- 
ments covering the circumference. Each of the segments 
was 7.5 mm thick and tt/6 radians wide between the 
endo- and epicardial border. Within each segment the 
material points were transformed to local Cartesian car- 
diac coordinates with X c defined as tangential to the 
surface of the volume and parallel with the short axis 
planes, X L defined as orthogonal to the short axis planes 
and oriented apex-to-base, and X R defined parallel with 
the short axis planes, orthogonal to X c and oriented 
outwards from the volume. 

Four orders of the polynomial function g(X) were ana- 
lyzed. 1) A bilinear-quadratic polynomial (BLQ) was 
bilinear within the X C -X L plane and quadratic in X R , 2) 
a bilinear-cubic polynomial (BLC) was bilinear within 
the X C -X L plane and cubic in X R , 3) a linear-quadratic- 
cubic polynomial (LQC) was cubic in X R , quadratic in 
X c and linear in X L and 4) a biquadratic-cubic polyno- 
mial (BQC) was biquadratic within the X C -X L plane and 
cubic in X R . For example, for the LQC polynomial field, 
the estimated coordinates were given by 

Xi = (auXl + a 2 iXl + a 3i X R + a Ai ) (a 5i X^ + a 6i X c + a 7i ) (a 8i X L + a 3i ) 



[* 3 rX 2 C Xl 1] 



Cli 



C24i 



(3) 



where a ki and c# are sets of constants, i = R, C, L, and 
k and / are serial numbers for the indices. The polyno- 
mial field for each segment was determined by minimiz- 
ing the squared difference between the measured and 
the estimated coordinates of the material points belong- 
ing to the segmentation. 



mm 



J2 ( Xi i 



(c,)) 2 , 



(4) 



where m is the total number of points within the seg- 
ment and Ci are all the constants c ti . The polynomial 
was evaluated and the Lagrangian strain tensor quanti- 
fied at material positions along the radial centerline 
within each segment. The accuracy of the polynomial 
function was evaluated by the residual of the polynomial 
mapping, determined as the root-mean-squared differ- 
ence (RMS) between the measured (x) and estimated 
(x) coordinates within each segment. The strain analysis 
was performed off line using Matlab 7.6.0 (The Math- 
Works, Inc.). 
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Analytical evaluation 

A previously presented analytical model [6,10] was used 
to evaluate the strain computation. Briefly, the model 
was a thick-walled cylinder that deforms to resemble the 
LV during systole. The cylinder was adapted to the sys- 
tolic in vivo data with inner radius of 28.5 mm and 
outer radius of 42.4 mm at the reference configuration, 
and inner radius of 23.6 mm at the deformed configura- 
tion. Material points were sampled through the cylinder 
wall within ten short axis planes with 2.5 mm separation 
between the planes and 2.5 x 2.5 mm separation 
between sampled points within each plane. 

Strains were computed throughout the model using 
the method described above. For all polynomial orders, 
the estimated strains were compared with the analytical 
strains (Mjj , where 7,/ = R,C,L). To mimic the acquisition 
of in vivo data, noise was added to the analytical model 
to investigate the noise sensitivity of the analysis 
method. 

Three normalized distributions of noise were analyzed, 
each with a mean value of zero and with the standard 
deviations 0.05 mm, 0.10 mm and 0.30 mm, respectively. 
The noise levels evaluated were chosen to represent 
practical measurement errors in DENSE displacement 
measurements corresponding to SNR in the order of 35, 
17 and 6. 

In vivo evaluation 

In vivo evaluation was performed by applying the strain 
quantification method to displacements acquired by 
DENSE MRI in a normal human heart, and comparing 
the results with strains obtained for normal human 
hearts from other studies. 

Anatomical images and displacement data were 
acquired in a 31-year old healthy volunteer on a 1.5 T 
MR scanner (Philips Achieva, Philips Medical Systems, 
Best, The Netherlands). Displacement data were 
acquired, using an in-house DENSE implementation, in 
three cardiac phases as illustrated in Figure 1; one in 
systole, and two during diastole related to early and late 
filling. The cardiac phases were defined from balanced 
steady state free precession acquisition in a standard 3 
chamber view with a temporal resolution of 13 ms. 
These images were used to define mitral valve opening 
(MVO) and end systole (ES), defined as aortic valve clo- 
sure (AVC). These two events were used to time the 
measurement of displacement in systole and diastole, 
acquired in 3D short axis slab in two separate 3D 
DENSE acquisitions. For the systolic measurement, the 
initial reference position was encoded at end diastole 
(ED), defined by the ECG R-wave peak, and the displa- 
cement was read out at ES, as shown in Figure 1. For 
the displacement in diastole, the initial reference posi- 
tion was encoded at MVO and read out at both MVO + 




Diastolic DENSE 



Figure 1 Timing of acquisition. Definitions of the time frames of 

end diastole (ED), end systole (ES), mitral valve opening (MVO) and 

the two time frames during diastolic filling (75 ms and 213 ms after 

MVO, respectively). ED is defined at the ECG R-wave peak and ES at 

aortic valve closure (AVC). 
k J 

75 ms and at MVO + 213 ms, which corresponds to 
approximately 22% and 62%, respectively, of the LV fill- 
ing interval. The 3D DENSE acquisitions were per- 
formed with an ECG triggered, respiratory navigator 
gated acquisition with the following parameters: acquisi- 
tion voxel size 2.5 x 2.5 x 2.5 mm, field of view 350 x 
350 x 25 mm, /c-space segmentation factor 3, EPI factor 
7, parallel imaging (SENSE) with reduction factor 2, 
echo time 5.7 ms, repetition time 12.2 ms, balanced 
multi-point displacement encoding 0.35 cycles per pixel 
[11], 3-SPAMM [12] displacement encoding with opti- 
mized flip angle [13]. The sequence requires data acqui- 
sition of 346 heart beats. Including navigator gating, the 
scan time is in the order of 10-15 minutes, depending 
on heart rate and navigator efficiency. For this DENSE 
sequence, a noise level of about 0.05 mm can be 
expected for a typical in-vivo measurement. 

Segmentation of the myocardium was performed using 
the freely available software Segment [14]. To avoid that 
non-myocardial tissue could affect the myocardial strain 
estimates, the segmentation was conservative, excluding 
uncertain areas or voxels. The phase of the DENSE MRI 
is proportional to the myocardial displacement, but can 
be an arithmetic modulo 2tt, which can be represented 
by wrapping the phase into the interval -tt to tt. Using 
the segmentation to outline the myocardium, the MRI 
phase was unwrapped [15] using a multigrid solver to 
the Poisson equations [16]. 

Systolic Lagrangian strains were analyzed at ES, with 
reference configuration at ED. Diastolic Lagrangian 
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strains were analyzed during LV filling with reference 
configuration at MVO. 

The research has been performed with informed con- 
sent, approved by the Regional Ethical Review Board in 
Linkoping and carried out in compliance with the Hel- 
sinki Declaration. 

Results 

Analytical evaluation 

The size of the errors of estimated strains is dependent 
on the spatial resolution of the sampled displacements 
[6]. Given the analytical strain components Ejj of the 
model, the absolute errors of the estimated strains in 
the analytical model, £ /7 = \^ ir Eu\ } were computed for 
each polynomial configuration order, integrated 
throughout each plane and averaged over the planes. 
For the noise free case, the mean ± SD error for the 
strain components of each polynomial order were BLQ: 
0.0098 ± 0.0067, BLC: 0.0125 ± 0.0092, LQC: 0.0068 ± 
0.0016, BQC: 0.0070 ± 0.0017. The errors for each com- 
ponent, respectively, from the linear-quadratic-cubic 
polynomial, for which the smallest mean errors were 
obtained, were s RR = 0.0081 for the radial strain, s cc = 
0.0054, s LL = 0.0089, s RC = 0.0072, s RL = 0.0048 and s CL 
= 0.0067. The LQC RMS differences, i.e. the residual of 
the polynomial mapping, averaged within the analytical 
model, were RMSx R = 0.030 mm in the radial direction, 
RMSx c = 0.017 mm in the circumferential direction, 
and RMSx L = 0.0004 mm in the longitudinal direction. 
The sensitivities to noise for each polynomial are plotted 
in Figure 2. The LQC polynomial yielded the smallest 
mean errors in the ideal situation and for a simulated 
noise with a standard deviation of 0.05 mm while the 
BLQ polynomial yielded the smallest mean error for the 
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Figure 2 Sensitivities to noise. The mean ± SD absolute error in 
estimated strain of each polynomial for different extents of noise 
on the analytical model. BLQ: bilinear-quadratic polynomial, BLC: 
bilinear-cubic polynomial, LQC: linear-quadratic-cubic polynomial, 
BQC: biquadratic-cubic polynomial. 



cases with 0.10 and 0.30 mm standard deviation of 
noise. 

In vivo evaluation 

The LQC and the BQC polynomials yielded the smallest 
maximum RMS differences in both systole and diastole. 
The results from the LQC polynomial are analyzed in 
further detail below. A mid-section (comprised of the 
short axis slices 5-7) is shown in Figure 3, 4, 5 and 6. 

Systole 

The systolic strain components estimated using the 
polynomial method are shown in Figure 3. For compari- 
son, Figure 4 shows the corresponding components esti- 
mated using a finite element method [17]. 

Maximum lengthening (up to 0.74) was found in the 
subendocardium in the septum and lateral free wall. 
Maximum shortening (down to -0.35) was found in the 
subendocardium and was essentially evenly distributed 
throughout the plane. 

For all three directions, the lowest RMS values were 
found in septum and the highest values in the posterior 
wall. RMSx R was within the range 0.13 - 0.40 mm, 
RMSx c 0.09 - 0.21 mm, and RMSx L 0.09 - 0.24 mm. 

Diastole 

All diastolic strain components at 75 ms (22% of the fill- 
ing interval) and 213 ms (62%) after MVO at the mid- 
section of the 3D volume are shown in Figure 5. 

The diastolic principal strains at 213 ms after MVO at 
the mid-section of the 3D volume are shown in Figure 
6. Maximum lengthening (up to 0.35) was found in the 
postero-lateral wall while maximum shortening was 
essentially evenly distributed over the plane. 

The highest RMS values were found in the posterior 
wall at 75 ms after MVO and in the anterior and poster- 
ior walls at 213 ms after MVO. The RMS values ranged 
from 0.07 mm to 0.40 mm for all directions at both dia- 
stolic times, except RMSx R at 75 ms after MVO which 
approached 0.50 mm in a small region in the posterior 
wall. 

Discussion 

The proposed myocardial strain quantification method 
was initially developed for strain computation on data 
from a surgically implanted transmural bead array. 
However, this work shows that the method can be 
extended to be used with displacement data from displa- 
cement-encoded MRI. 

This work uses a polynomial function to find a differ- 
entiable expression from the discrete displacement field. 
This polynomial function assumes continuous displace- 
ment in the myocardium, which reflects the connective 
properties of the myocardial tissue. This a priori 
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Erc Erl Ecl 

Figure 3 Systolic strains. All components of the end-systolic strain tensor at the mid-section of the 3D volume in a healthy volunteer 
estimated using the polynomial method. 




Erc Erl Ecl 

Figure 4 Systolic strains. All components of the end-systolic strain tensor at the mid-section of the 3D volume in a healthy volunteer 
estimated using a finite element method. 
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Figure 5 Diastolic strains. All components of the diastolic Lagrangian strain tensor at the mid-section of the 3D volume in a healthy volunteer, 
a) At 75 ms after MVO; b) At 213 ms after MVO. 



information helps making the estimation more robust to along each dimension, which in turn depends on the 

noise. spatial resolution of the data, wall thickness and the 

The optimal order of the polynomial functions in sizes of the finite segments. Generally, the number of 

equation (1) depends on the number of material points unknown constants in the polynomial should be less 
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Figure 6 Diastolic principal strains. Diastolic principal strains at the mid-section of the 3D volume in a healthy volunteer, at 213 ms after MVO. 
a) E1; b) E2; c) E3. 



than the number of measured points along each dimen- 
sion, which implies that a third order polynomial 
requires at least five measured material points along the 
corresponding dimension, a second order polynomial 
requires at least four points and a first order polynomial 
requires at least three points. Furthermore, to avoid an 
undetermined problem, the number of data points must 
be equal or greater than the number of coefficients in 
the minimization. This implies that the minimum num- 
ber of required data points depends on the polynomial 
orders; BLQ requires 12 data point, BLC 16, LQC 24, 
and BQC 36. Hence, the polynomial order is limited by 
the number of included data points. This requirement 
was fulfilled for the in vivo evaluation; however, the 
margin for the BQC polynomial was small at some loca- 
tions of the myocardium. 

Four different polynomial orders were evaluated. The 
smallest absolute errors of the estimated strains in the 
analytical model in the presence of low noise were 
obtained with a linear-quadratic-cubic polynomial. In 
the subsequent in vivo validation, the in vivo results of 
the linear-quadratic-cubic polynomial were analyzed in 
detail. Given the incumbent spatial resolution, the 
restriction on fR(X R ) was the wall thickness at the late 
diastolic time frame and the restrictions on fc(X c ) and 
f L (X L ) were the width and height of each segment, 
respectively. The width (tt/6 radians) and height (7.5 
mm) of each segment were kept small in order to 
resolve local variations of deformation. 

The RMS differences between the acquired in vivo 
coordinates and the coordinates estimated using the 
polynomial fitting method reflect the accuracy of the 
polynomial fitting, giving a comprehension of the extent 
to which the coordinates estimated by the polynomial fit 
to the acquired in vivo coordinates. For both systolic 
and diastolic data, the RMS differences of the LQC 
polynomial were highest in the posterior wall. For the 
systolic data the region of highest RMS differences was 



close to the posterior papillary muscle, which might 
have caused locally irregular displacements. The RMS 
differences may furthermore be affected by the spatially 
varying signal-to-noise ratio in the DENSE measurement 
[18]. For the diastolic data, the region of highest RMS 
differences coincided with the region of thinnest wall. 
The present implementation of the method used the 
same polynomial order for each segment within the 3D 
volume. If a data set with large variations in wall thick- 
ness is considered, e.g. a 3D slab comprising an infarct 
with thin wall near the infarcted myocardium, the high- 
est accuracy could be obtained by local adjustment of 
the size of the segments to the size of the infarcted 
region and adapting the polynomial order to wall thick- 
ness and segment width. 

Systolic radial, circumferential and longitudinal strains, 
as well as systolic circumferential-longitudinal shear, 
show agreement with systolic strains previously reported 
for human myocardium [19,20]. E RC : We observed 
somewhat higher magnitudes of radial-circumferential 
shear strain than the results of Moore et al. [19], and we 
observed the lowest values in the anterior region of the 
myocardium while Moore et al. reported the lowest 
values in septum. E RL : The observed radial-longitudinal 
shear values fits within mean ± 2SD of the values 
reported by Moore et al. [19]. 

Diastolic function of the LV is determined by a com- 
plex sequence of many interrelated events and para- 
meters including active relaxation, elastic recoil, passive 
filling characteristics, heart rate and inotropic state. Dia- 
stolic LV filling is a highly dynamic process with early 
and late transmural inflows. Thus a detailed analysis of 
myocardial strain during diastole requires resolving the 
temporal process of diastolic filling. 

The highest values of circumferential strain during the 
first 213 ms of diastolic filling were observed in the pos- 
tero-lateral wall. The same quantitative behavior has 
been reported in previous studies of early diastolic 
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strains in normal human hearts [21,22]. Radial strain, 
interpreted as wall thinning during diastole, was most 
apparent in the lateral wall, which also has been 
reported by others [21]. 

Limitations 

This work is limited to study the kinematics of the 
heart, focusing on strain. Strain should preferably be 
related to an unloaded, stress free reference configura- 
tion. Using in vivo data, there is no unloaded, stress free 
configuration of the heart. Instead, the reference config- 
urations correspond to defined time points in the car- 
diac cycle. The strain presented in this, and similar 
articles within the field, therefore disregards the residual 
strain needed to study cardiac kinetics as opposed to 
kinematics. 

In vivo validation of strain is challenging, which is 
why an analytical model was included in the valida- 
tion. The analytical model can however never fully 
describe the cardiac kinematics. For in vivo estima- 
tion, the quality of the strain measurements is highly 
dependent on the quality of the underlying displace- 
ment data. While the polynomial fit reduces sensitiv- 
ity-to-noise in the measurements, image artifacts or 
errors in the displacement measurements can deterio- 
rate the strain estimates. Improving the DENSE acqui- 
sitions is an active field in the MRI research 
community, and strain analysis methods like the one 
presented in this paper will benefit directly from such 
improvements. 

Conclusions 

In conclusion, the proposed method for strain estima- 
tion within a 3D myocardial volume yields accurate 
results when validated on an analytical model. The poly- 
nomial field is capable of resolving the measured mate- 
rial positions from the in vivo data, and the obtained in 
vivo strains agree with previously reported myocardial 
strains in normal human hearts. 
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